import os, math
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.offline as pyo
import matplotlib.pyplot as plt
%matplotlib inline
import plotly.io as pio
pio.renderers.default = "notebook+pdf"
data_dict = {}
header = {}
i = 0
for filename in os.listdir():
if "_data" in filename:
if "_spect" in filename:
data_dict[filename[:-15]] = pd.read_csv(filename)
else:
data_dict[filename[:-9]] = pd.read_csv(filename)
elif "_header" in filename:
head = pd.read_csv(filename, names = ["Name", "Value"])
if "_spect" in filename:
header[filename[:-17]] = head.at[16,"Value"]
else:
header[filename[:-11]] = head.at[16,"Value"]
i += 1
print("Example of dataframe:")
data_dict["background"]
Example of dataframe:
| Channel | Energy | Counts | |
|---|---|---|---|
| 0 | 0 | -62.744125 | 0 |
| 1 | 1 | -60.751958 | 0 |
| 2 | 2 | -58.759791 | 0 |
| 3 | 3 | -56.767624 | 0 |
| 4 | 4 | -54.775457 | 0 |
| ... | ... | ... | ... |
| 1019 | 1019 | 1967.274151 | 3 |
| 1020 | 1020 | 1969.266319 | 1 |
| 1021 | 1021 | 1971.258486 | 4 |
| 1022 | 1022 | 1973.250653 | 1 |
| 1023 | 1023 | 1975.242820 | 1 |
1024 rows × 3 columns
import plotly.graph_objects as go
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict["180s_2017_sample6"]["Energy"],y=data_dict["180s_2017_sample6"]["Counts"],
mode='lines',
name="180s_2017_sample6"))
fig3.add_trace(go.Scatter(x=data_dict["BlockT_180s_2017_sample6"]["Energy"],y=data_dict["BlockT_180s_2017_sample6"]["Counts"],
mode='lines',
name='BlockT_180s_2017_sample6'))
fig3.add_trace(go.Scatter(x=data_dict["background"]["Energy"],y=data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
mode='lines',
name='Background'))
#fig3.update_yaxes(type="log")
def ECalib(ener):
return (ener+62.74412532637075)/1.9921671018276763
def integPeak(xmin,xmax,df):
minIn = round(ECalib(xmin))
maxIn = round(ECalib(xmax))
return df["Counts"][minIn:maxIn].sum()
def halflife(N0, Nt, t):
return t/math.log(Nt/N0,0.5)
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict["2013_sample5"]["Energy"],y=data_dict["2013_sample5"]["Counts"]*float(header["2017_sample7"])/float(header["2013_sample5"]),
mode='lines',
name='2013_sample5'))
fig3.add_trace(go.Scatter(x=data_dict["2017_sample7"]["Energy"],y=data_dict["2017_sample7"]["Counts"],
mode='lines',
name='2017_sample7'))
fig3.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
p1 = integPeak(1120.603,1246.11,data_dict["2013_sample5"]*float(header["2017_sample7"])/float(header["2013_sample5"]))
p1p = integPeak(1120.603,1246.11,data_dict["2017_sample7"])
print("Integral of peak 1: "+str(p1))
print("Integral of peak 1': "+str(p1p))
p2 = integPeak(1256.07,1437.58,data_dict["2013_sample5"]*float(header["2017_sample7"])/float(header["2013_sample5"]))
p2p = integPeak(1256.07,1437.58,data_dict["2017_sample7"])
print("Integral of peak 2: "+str(p2))
print("Integral of peak 2': "+str(p2p))
print()
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
print()
print("HL from peak (-1): " + str(halflife(p2p, p2, 3)))
print("HL from peak (0): " + str(halflife(p2p, p2, 4)))
print("HL from peak (+1): " + str(halflife(p2p, p2, 5)))
Integral of peak 1: 16623.9933101173 Integral of peak 1': 30989 Integral of peak 2: 13022.771535923754 Integral of peak 2': 23920 HL from peak (-1): 3.338938216649393 HL from peak (0): 4.451917622199191 HL from peak (+1): 5.564897027748988 HL from peak (-1): 3.4200472219342024 HL from peak (0): 4.560062962578936 HL from peak (+1): 5.700078703223671
s2013 = "2013_sample6"
s2017 = "2017_sample8"
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*float(header[s2017])/float(header[s2013]),
mode='lines',
name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
mode='lines',
name=s2017))
fig3.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
#title='Measure of {} and {} Gamma Ray Production, Counts vs. Energy'.format(s2013,s2017),
xmin = 765.9974
xmax = 925.3708
p1 = integPeak(xmin,xmax,data_dict["2013_sample6"]*float(header["2017_sample8"])/float(header["2013_sample6"]))
p1p = integPeak(xmin,xmax,data_dict["2017_sample8"])
print("Integral of peak 1: "+str(p1))
print("Integral of peak 1': "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
Integral of peak 1: 1711.284378109453 Integral of peak 1': 9781 HL from peak (-1): 1.1928892120792631 HL from peak (0): 1.5905189494390177 HL from peak (+1): 1.988148686798772
print("HL from p1: " + str(halflife(276, 35.12, 4)))
HL from p1: 1.3448526611739258
s2013 = "2013_sample7"
s2017 = "2017_sample1"
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
mode='lines',
name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
mode='lines',
name=s2017))
fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
fig3.show()
######################## Taking the integrals of the peak and calculating Half-Life #########################
xmin = 222.1358
xmax = 423.3446
p1 = integPeak(xmin,xmax,data_dict[s2013]*(float(header[s2017])/float(header[s2013])))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
First integral: 197111.41901922962 Second integral: 299631 HL from peak (-1): 4.965444340086574 HL from peak (0): 6.620592453448765 HL from peak (+1): 8.275740566810956
#Half life from single point
print("HL from p1: " + str(halflife(7604, 4958.303, 4)))
HL from p1: 6.4839059410690885
s2013 = "2013_sample8"
s2017 = "2017_sample5"
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
mode='lines',
name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
mode='lines',
name=s2017))
#fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
fig3.show()
######################## Taking the integrals of the peak and calculating Half-Life #########################
xmin = 1054.862
xmax = 1162.439
p1 = integPeak(xmin,xmax,data_dict[s2013]*float(header[s2017])/float(header[s2013]))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
First integral: 516.3999696140991 Second integral: 2055 HL from peak (-1): 1.5055875568362662 HL from peak (0): 2.007450075781688 HL from peak (+1): 2.50931259472711
#Half life from single point
print("HL from p1: " + str(halflife(62, 6.54, 4)))
HL from p1: 1.232701471488777
s2013 = "2013_sample4"
s2017 = "2017_sample4"
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict[s2013]["Energy"],y=data_dict[s2013]["Counts"]*(float(header[s2017])/float(header[s2013])),
mode='lines',
name=s2013))
fig3.add_trace(go.Scatter(x=data_dict[s2017]["Energy"],y=data_dict[s2017]["Counts"],
mode='lines',
name=s2017))
fig3.update_yaxes(type="log")
fig3.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
fig3.show()
######################## Taking the integrals of the peak and calculating Half-Life #########################
xmin = 56.7859
xmax = 118.5431
p1 = integPeak(xmin,xmax,data_dict[s2013]*float(header[s2017])/float(header[s2013]))
p1p = integPeak(xmin,xmax,data_dict[s2017])
print("First integral: "+str(p1))
print("Second integral: "+str(p1p)+"\n")
print("HL from peak (-1): " + str(halflife(p1p, p1, 3)))
print("HL from peak (0): " + str(halflife(p1p, p1, 4)))
print("HL from peak (+1): " + str(halflife(p1p, p1, 5)))
First integral: 6595.638372747748 Second integral: 54951 HL from peak (-1): 0.9808532444341904 HL from peak (0): 1.3078043259122538 HL from peak (+1): 1.6347554073903172
shields_Al = ["G","J","M","O"]
integrAl = {}
thicknAl = {"G":0.508,"J":1.016,"M":2.032,"O":2.54}
surfDensAl = {"G":141,"J":258,"M":522,"O":655}
fig = go.Figure()
for shield in shields_Al:
fig.add_trace(go.Scatter(x=data_dict["Block"+shield+"_130s_2017_sample3"]["Energy"],y=data_dict["Block"+shield+"_130s_2017_sample3"]["Counts"],
mode='lines',
name="Block"+shield))
integrAl[shield] = integPeak(580.7258,760.0209,data_dict["Block"+shield+"_130s_2017_sample3"])
fig.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
fig.show()
print(integrAl)
{'G': 9984, 'J': 10059, 'M': 9685, 'O': 9765}
shields_Pb = ["Q","R","S","T"]
integrPb = {}
thicknPb = {"Q":0.8128,"R":1.6256,"S":3.81,"T":6.35}
surfDensPb = {"Q":1120,"R":2066,"S":3448,"T":7367}
fig = go.Figure()
for shield in shields_Pb:
fig.add_trace(go.Scatter(x=data_dict["Block"+shield+"_130s_2017_sample3"]["Energy"],y=data_dict["Block"+shield+"_130s_2017_sample3"]["Counts"],
mode='lines',
name="Block"+shield))
integrPb[shield] = integPeak(580.7258,760.0209,data_dict["Block"+shield+"_130s_2017_sample3"])
fig.update_layout(xaxis_title='Energy [KeV]',
yaxis_title='Counts')
fig.show()
print(integrPb)
{'Q': 9144, 'R': 8521, 'S': 7131, 'T': 4923}
integr = {**integrAl,**integrPb}
thickn = {**thicknAl,**thicknPb}
surfDens = {**surfDensAl,**surfDensPb}
testArr = {'Surface Density [g/cm]': [], 'Thickness [mm]': [],'Counts in Peak': []}
a = []
for letter in shields_Al+shields_Pb:
testArr['Surface Density [g/cm]'].append(surfDens[letter])
testArr['Thickness [mm]'].append(thickn[letter])
testArr['Counts in Peak'].append(integr[letter])
shieldDf = pd.DataFrame(testArr)
shieldDf
| Surface Density [g/cm] | Thickness [mm] | Counts in Peak | |
|---|---|---|---|
| 0 | 141 | 0.5080 | 9984 |
| 1 | 258 | 1.0160 | 10059 |
| 2 | 522 | 2.0320 | 9685 |
| 3 | 655 | 2.5400 | 9765 |
| 4 | 1120 | 0.8128 | 9144 |
| 5 | 2066 | 1.6256 | 8521 |
| 6 | 3448 | 3.8100 | 7131 |
| 7 | 7367 | 6.3500 | 4923 |
from plotly.subplots import make_subplots
import seaborn as sns
fig3 = make_subplots(rows=1, cols=2, subplot_titles=("Default Scale", "Logarithmic Scale"))
fig3.add_trace(go.Scatter(x=shieldDf['Surface Density [g/cm]'],y=shieldDf['Counts in Peak'],
mode='markers',
name="180s_2017_sample6"),col=1,row=1)
fig3.add_trace(go.Scatter(x=shieldDf['Thickness [mm]'],y=shieldDf['Counts in Peak'],
mode='markers',
name='BlockT_180s_2017_sample6'),col=2,row=1)
fig = px.scatter(shieldDf, x='Surface Density [g/cm]', y='Counts in Peak', trendline="ols")
fig.show()
results = px.get_trendline_results(fig)
results.px_fit_results.iloc[0].summary()
/home/jonas/anaconda3/lib/python3.7/site-packages/scipy/stats/stats.py:1535: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=8
| Dep. Variable: | y | R-squared: | 0.986 |
|---|---|---|---|
| Model: | OLS | Adj. R-squared: | 0.984 |
| Method: | Least Squares | F-statistic: | 436.1 |
| Date: | Wed, 18 Nov 2020 | Prob (F-statistic): | 7.85e-07 |
| Time: | 00:08:52 | Log-Likelihood: | -53.558 |
| No. Observations: | 8 | AIC: | 111.1 |
| Df Residuals: | 6 | BIC: | 111.3 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| const | 1.007e+04 | 104.686 | 96.152 | 0.000 | 9809.595 | 1.03e+04 |
| x1 | -0.7263 | 0.035 | -20.883 | 0.000 | -0.811 | -0.641 |
| Omnibus: | 4.563 | Durbin-Watson: | 2.390 |
|---|---|---|---|
| Prob(Omnibus): | 0.102 | Jarque-Bera (JB): | 1.403 |
| Skew: | -1.016 | Prob(JB): | 0.496 |
| Kurtosis: | 3.277 | Cond. No. | 3.95e+03 |
from plotly.subplots import make_subplots
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=data_dict["180s_2017_sample6"]["Energy"],y=data_dict["180s_2017_sample6"]["Counts"]-data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
mode='lines',
name="180s_2017_sample6"))
fig3.add_trace(go.Scatter(x=data_dict["BlockT_180s_2017_sample6"]["Energy"],y=data_dict["BlockT_180s_2017_sample6"]["Counts"]-data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
mode='lines',
name='BlockT_180s_2017_sample6'))
fig3.add_trace(go.Scatter(x=data_dict["background"]["Energy"],y=data_dict["background"]["Counts"]*float(header["180s_2017_sample6"])/float(header["background"]),
mode='lines',
name='Background'))
fig3.update_yaxes(type="log")